#ifndef POINT_H
#define POINT_H
#include <math.h>
#include <string.h>
#include <iomanip>
#include <iostream>
template <typename real>
class Latt;
template <typename real>
class point
{
private:
    real P[3];
public:
    point(){
        memset(P,0,sizeof(real)*3);
    }
    point(const real &orig){
        memcpy(P,orig.P,sizeof(real)*3);
    }
    point(const real &a,const real &b, const real &c){
        P[0]=a;P[1]=b;P[2]=c;
    };
    /*void operator=(const point &orig){
        memcpy(P,orig.P,sizeof(real)*3);
    }*/
    point operator- () const{
        point tmp;
        tmp.P[0]=-P[0];
        tmp.P[1]=-P[1];
        tmp.P[2]=-P[2];
        return tmp;
    }
    real &operator[](const int &i){
        return P[i];
    }
    real  operator[](const int &i)const {
        return P[i];
    }
    point operator+ (const point &rhs) const {
        point tmp; 
        tmp.P[0]=P[0]+rhs.P[0];
        tmp.P[1]=P[1]+rhs.P[1];
        tmp.P[2]=P[2]+rhs.P[2];
        return tmp;
    }
    point operator- (const point &rhs) const{
        point tmp;
        tmp.P[0]=P[0]-rhs.P[0];
        tmp.P[1]=P[1]-rhs.P[1];
        tmp.P[2]=P[2]-rhs.P[2];
        return tmp;
    }
    real  operator* (const point &rhs) const {
        return P[0]*rhs.P[0]+P[1]*rhs.P[1]+P[2]*rhs.P[2];
    }
    point operator* (const real   rhs) const {
        point tmp;
        tmp.P[0]=P[0]*rhs;
        tmp.P[1]=P[1]*rhs;
        tmp.P[2]=P[2]*rhs;
        return tmp;
    }
    point operator/ (const real   rhs) const {
        point tmp;
        tmp.P[0]=P[0]/rhs;
        tmp.P[1]=P[1]/rhs;
        tmp.P[2]=P[2]/rhs;
        return tmp;
    }
    void  operator*=(const real  time){
        P[0]*=time;
        P[1]*=time;
        P[2]*=time;
    }
    void  operator/=(const real   rhs){
        P[0]/=rhs;
        P[1]/=rhs;
        P[2]/=rhs;
    }
    void  operator+=(const point &orig){
        P[0]+=orig.P[0];
        P[1]+=orig.P[1];
        P[2]+=orig.P[2];
    }
    void  operator-=(const point &orig){
        P[0]-=orig.P[0];
        P[1]-=orig.P[1];
        P[2]-=orig.P[2];
    }   
    void  operator*=(const Latt<real> &);
    real  norm() const {
        return sqrt(P[0]*P[0]+P[1]*P[1]+P[2]*P[2]);
    }
    void show(int n=6){
        std::cout << std::fixed << std::setprecision(n)<<std::setw(n+4) <<P[0]<<std::setw(n+5)<<P[1]<<std::setw(n+5)<<P[2];
    }
};

template <typename real>
point<real> operator*(const real &lhs,const point<real> &rhs)
{
    point<real> tmp;
    tmp[0]=lhs*rhs[0];
    tmp[1]=lhs*rhs[1];
    tmp[2]=lhs*rhs[2];
    return tmp;
}

template <typename real>
point<real> cross(const point<real> &lhs,const point<real> &rhs)
{
    point<real> tmp;
    tmp[0]=lhs[1]*rhs[2]-lhs[2]*rhs[1];
    tmp[1]=lhs[2]*rhs[0]-lhs[0]*rhs[2];
    tmp[2]=lhs[0]*rhs[1]-lhs[1]*rhs[0];
    return tmp;
}

template <typename real>
real cos(const point<real> &lhs,const point<real> &rhs)
{
    return (lhs*rhs)/(lhs.norm()*rhs.norm());
}
#endif